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ABSTRACT 

•/ 

The  resistance  and  mobility  functions  which  completely  characterize  the 
linear  relations  between  the  force,  torque  and  stresslet  and  the  translational 
and  rotational  velocities  of  two  spheres  in  low-Reynolds-number  flow  have  been 
calculated  using  a  boundary  collocation  technique.  The  ambient  velocity  field 
is  assumed  to  be  a  superposition  of  a  uniform  stream  and  a  linear  (vorticity 
and  rate-of-strain)  field.  This  is  the  first  compilation  of  accurate 
expressions  for  the  entire  set  of  functions.  Our  calculations  are  in 
agreement  with  earlier  results  for  all  functions  for  which  such  results  are 
available.  Our  technique  is  successful  at  all  sphere-sphere  separations 
except  at  the  almost-touching  (gaps  of  less  than  .005  diameter)  configuration. 

<L 

New  results  for  the  stresslet  functions  have  been  used  to  determine 
Batchelor  and  Green's  (1972)  order  .C2  coefficient  in  the  bulk-stress  (7.1 
instead  of  their  7.6).  The  two-sphere  functions  have  also  been  used  to 
determine  the  motion  of  a  rigid  dumbbell  in  a  linear  field.  We  also  show  that 
certain  functions  have  extrema.  The  source  (FORTRAN)  code  is  furnished  in  the 
appendix,  j 
/ 

/ 
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SIGNIFICANCE  AND  EXPLANATION 


v: 


The  calculation  of  hydrodynamic  interactions  between  particles  ia  needed 
for  the  understanding  and  control  of  many  natural  and  manufacturing  processes, 
for  instance,  those  involving  sedimentation,  colloidal  stability,  ^aspension 
rheology,  and  cloud  formation.  A  fundamental  approach  to  these  problems  often 
requires  detail  information  on  the  hydrodynamic  interactions  between  two 
spheres,  that  is,  the  forces,  torques  and  stress  dipoles  induced  by  the 
particle  motions  and  the  ambient  velocity.  Until  now,  the  available 
information  was  incomplete. 

This  report  furnishes  the  complete  solution  of  the  problem  using  a 
collocation  approach.  The  results  are  in  excellent  agreement  with  all  earlier 
solutions  for  special  cases  of  the  complete  problem.  The  source  (FORTRAN) 
code  has  been  included. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC  and  not  with  the  authors  of  this  report. 
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4.  Results  and  Conclusions 


-  results  for  the  resistance  and  mobility  problems 

-  comparisons  in  the  far-field  with  the  method  of  reflections 

-  comparison  with  the  method  of  twin-multi pole  expansions 

-  comparison  with  solutions  by  bispherical  co-ordinates 

-  comparisons  in  the  near-field  with  lubrication  solutions 


5.  Sample  Problems 

2 

5.1  The  Bulk  Stress  in  a  Suspension  of  Spheres  to  Order  c 

5.2  Motion  of  a  Rigid  Dumbbell  in  a  Shear-Field 


Appendix  1 .  Notes  on  the  Computer  Programs 


Appendix  2.  Program  Listings 
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THE  RESISTANCE  AND  MOBILITY  FUNCTIONS  OF  TWO  EQUAL 

SPHERES  IN  LOW  REYNOLDS  NUMBER  FLOW 

*  1  2  **13 

Sangtae  Kim  *  ’  and  Richard  T.  Mifflin  *  * 

1 .  Introduction 

In  this  work  we  have  used  the  boundary  collocation  technique  to  calculate  the 
set  of  functions  which  describe  the  hydrodynamic  interaction  between  two  rigid 
spheres  in  low-Reynolds-number  flow.  Such  information  is  needed  in 
theoretical  investigations  of  the  behavior  of  suspensions  of  small 
(sub-micron)  particles  as  shown  in  the  review  articles  by  Batchelor  ( 1 97^)  and 
Jeffrey  and  Acrivos  (1976).  Specific  applications  are  found  in  studies  of 
sedimentation  velocities  (Batchelor  1972),  rheological  properties  (Batchelor 
and  Green  1972  ),  Brownian  diffusion  (Batchelor  1976)  and  fixed-bed 
permeabilities  (Howells  1974).  In  all  cases  the  specific  Information  that  is 
required  is  the  linear  relation  between  the  rigid-body  motion  of  the  spheres, 

ya  +  80x(-’?a)’ 

in  an  ambient  field,  U(x)  -  U*  ♦  Q*x  *  E*x  on  one  hand  and  on  the  other  hand 

the  force,  torque  and  stresslet  (moment)  exerted  by  each  sphere  on  the  fluid. 

The  notation  is  as  follows:  U  and  0  are  the  translational  and  rotational 

-a  -o 

velocities  of  the  sphere  centered  at  x  ,  o-1,2;  u",  Q  and  E  are  the  uniform 

-a  -  -  _ 

stream,  constant  vorticity  and  rate-of-straln  fields.  The  force,  torque  and 
stresslet  on  each  sphere  are  given  by  the  following  integrals  of  the  stress, 
o,  over  the  surface  of  sphere  a: 


F 

-a 


-/c  o*n  dA, 
P 


T  -  —  f q  (x-x  )x(o*n)  dA, 
-a  -  -0  ■ 
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-a  "  2^S  )o*n  +  2*n(x-x  )]dA. 

P  U  m  —  — 


n  is  the  outward  normal  vector  for  the  surface. 

Because  of  the  linearity  of  the  governing  equations,  the  solution  of  our 
hydrodynamic  interaction  problem  (with  equal  spheres)  is  completely  specified 
by  19  independent  scalar  functions.  We  present  here  what  we  believe  to  be  the 
first  complete  solution  of  this  problem.  We  emphasize  that  our  method  can  be 
reduced  to  a  60-line  routine  that,  with  the  help  of  subprograms  for  the 
special  functions  (Legendre  functions)  calculates  the  entire  collection  of 
functions.  Others  have  solved  various  subsets  of  this  problem  using 
bispherical  coordinates  (Stimson  and  Jeffery  1926;  Goldman,  Cox  and  Brenner 
1966;  Lin,  Lee  and  Sather  1970),  method  of  reflections  (Happel  and  Brenner 
1965)  and  lubrication  theory  (O'Neill  and  Majumdar  1970),  as  reviewed  by 
Jeffrey  and  Onlshi  ( 1 98^0 .  These  authors  also  present  a  comprehensive 
solution  of  the  important  subset  involving  the  force  and  torque  in  an  ambient 
field  composed  of  a  uniform  stream  and  vorticity  field,  using  the 
twin-multipole  variation  of  the  method  of  reflections. 

A  more  detailed  discussion  of  hydrodynamic  interaction  is  presented  in 


following  subsections  where  we  review  the  resistance  and  mobility  functions. 
In  Section  2,  we  show  how  the  boundary  collocation  technique  of  Gluckman 


r 


et.  al.  (1971)  and  Lamb's  velocity  representation  can  be  used  to  solve  the 
boundary  value  problem  associated  with  each  resistance  function.  Formulae 
which  relate  the  resistance  functions  to  the  coefficients  in  the  velocity 
representation  are  found  in  Section  3.  We  state  our  principal  results  in 
Section  M  and  illustrate  sample  applications  in  Section  5,  including  the 
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correction  of  Batchelor  and  Green's  (1972b)  result  for  the  coefficient  of  the 
2 

c  term  in  the  bulk  stress  (7.1  instead  of  their  7.6)  and  the  hydrodynamic 
functions  for  the  rotation  of  a  rigid  dumbbell  in  a  linear  field.  We  have 
placed  our  source  code  and  sample  calculations  in  the  appendix. 
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1.1  The  Resistance  Problem 


Following  Brenner  and  O'Neill  (1972),  we  define  the  resistance  problem  as  that 
in  which  the  force,  torque  and  stresslet  are  to  be  determined  for  a  specified 
particle  motion  in  the  ambient  field.  The  linearity  of  the  Stokes  equations 
permits  the  expression  of  the  forces,  torques  and  stresslets  in  the  following 
matrix  form: 
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The  6X6  matrix  of  tensors  has  been  named  the  grand  resistance  matrix  by 
Rallison  (1977). 


We  have  followed  the  development  and  notation  in  Jeffrey  and  Onishi 
(198i|)  and  Jeffrey  ( 1 984)  throughout  this  section.  A,  B,  B  and  C  are  second 
rank  tensors,  G,  G,  H  and  H  are  third  rank  tensors  and  the  M  tensors  are 
fourth  rank  tensors.  We  shall  see  below  that  there  are  inter-relations 
between  certain  tensor  pairs.  These  pairings  are  highlighted  by  using  the 
same  letters  and  the  symbol,  -. 

We  first  reduce  the  number  of  independent  tensors  by  using  properties 
that  are  independent  of  the  particle  geometry.  The  reciprocal  theorem  of 
Lorentz  (1906)  can  be  used  to  show  that  the  resistance  matrix  is  symmetric 
(Brenner  and  O'Neill  1972  and  Hlnch  1972),  i.e.: 


(1 .1a,b,c) 
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We  may  impose  additional  relations  because  S  and  E  are  symmetric  and 
traceless.  The  condition  on  S  permit  us  to  set 
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on  E  require  that 
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The  symmetry  of  the  two-sphere  geometry  implies  that  each  tensor 
satisfies: 


Pa6(R)  -  P(3_o)(3“6)(-R)f 


(1.4) 


where  R=x2~x1  is  the  center- to-center  vector.  Finally,  the  axisymmetry  about 
the  sphere-sphere  axis  implies  that  each  tensor  can  be  decomposed  into 
expressions  involving  no  more  than  three  scalar  functions  (Brenner  1963, 

P 

1964).  Jeffrey  and  Onishi  (1984)  designate  these  scalar  functions  as  X  D(R), 

ap 

p  p 

Yag(R)  and  Z^tR)  with  P,  a  and  &  denoting  the  appropriate  tensor  P  .  They 
reserve  the  letter  X  for  those  funtions  that  arise  from  axisymmetric  flows. 
(More  specifically,  we  shall  see  later  that  the  X,  Y,  and  Z  functions  arise 
from  boundary  conditions  involving  spherical  harmonics  with  the  azimuthal 
number,  m,  equal  to  0,  1  and  2  respectively).  Thus,  with  d  -  R/R: 
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(1 . 5e) 


Ci.5f) 
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We  now  nondimensionalize  these  scalar  functions  so  that  they  become 
functions  only  of  the  dimensionless  separation  parameter  R/a.  The 


dimensionless  functions  will  be  denoted  with  the  symbol 


i“8  -  6,aAae.  BaS  -  k,a2ia6,  coB  -  8»a3CoB, 

GaB  -  «™2CaB.  HoB  -  8,a3KoB,  HoB  -  f^A"6. 


(1 .6a-c) 


(1.6d-f) 


The  dimensionless  functions  for  the  tensors  on  the  diagonal  of  the  grand 
resistance  matrix  will  approach  unity  for  large  R  because  the  scales  were 
chosen  by  considering  the  single-sphere  result. 


1 .2  The  Mobility  Problem 

Following  Batchelor  (1976),  we  define  mobility  problems  as  those  in  which  the 
particle  forces  and  torques  are  prescribed  in  the  ambient  field  and  the 
particle  motion  and  stresslet  are  the  unknowns.  The  formulation  of  the 
mobility  problem  is  rather  awkward  from  a  mathematical  perspective  since  the 
boundary  conditions  involve  the  unknowns,  but  in  many  problems  the  forces  and 
torques  are  the  prescribed  physical  quantities  and  the  particles  must  move 
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accordingly.  In  later  sections,  we  shall  first  solve  the  resistance  problem 
and  then  use  the  relations  between  the  mobility  and  resistance  tensors  to 


solve  the  mobility  problem.  Again,  the  linearity  of  the  Stokes  equation 
allows  us  to  write: 
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As  in  the  resistance  problem,  the  number  of  unknowns  can  be  reduced  by 
applying  the  reciprocal  theorem,  and  the  consequences  of  S  and  E  being 
symmetric  and  traceless.  Thus: 
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(1 .7a,b,c) 
(1.7d,e) 
(1.7f) 


and  equations  which  are  analagous  to  (1.2)  and  (1.3). 

As  before,  the  two-sphere  symmetry  allows  the  following  decompositions: 
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The  nondimensionalizatlons  of  these  functions  are  as  follows: 


,“8  -  6,aa“‘, 
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c“8  -  «,aV’. 

(1 .9a-c) 

-  Ba8/(2a). 
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(1 .9d-f ) 

1 .3  Relations  between  the  Resistance  and  Mobility  Functions 
Our  numerical  technique  solves  the  resistance  problem.  We  obtain  the  mobility 
functions  by  using  the  following  relations  between  the  (dimensional)  mobility 
and  resistance  functions.  In  matrix  form,  we  have: 


I 


.8 

11 

y8 

y12 

h 

*y11 

h 

"y12 

.6 

21 

yg 

y22 

h 

-y21 

h 

"y22 

r  yg 

ii 

yg 

12 

-yh 

X11 

Y° 

L  21 

yg 

x22 

-yh 

21 

va  a 

b 

b 

y11  -12 

y11 

y21 

a  a 

b 

b 

y21  y22 

y1 2 

y22 

b  b 

6 

c 

cvj 

>T 

r— 

>> 

y11 

y1 2 

fa 

„c 

c 

y21  y22 

y21 

y22 

♦  XG  (x8 
o2v  12 

♦  x8  ) 
*22' 

9 

+  YG  (y8 
xa2'y12 

+  y|2) 

♦  Y*  (yh 

xo2^y12 

9 

(1.14) 


(1.15a-c) 


"(Za1  +  Zo2)’  for  a  “  1’2* 


The  above  equations  hold  for  two  unequal  spheres.  If  we  limit  the 
analysis  to  the  case  of  two  equal  spheres,  then  the  symmetry  relation, 
equation  (1.4),  implies  that  subscripts  "22"  and  n21”  may  be  replaced 
everywhere  by  "11"  and  "12"  respectively.  For  functions  associated  with  B,  b, 
G  and  g,  this  substitution  requires  a  change  in  sign.  From  here  on,  without 
loss  of  generality,  we  shall  restrict  our  attention  to  the  "11"  and  "12" 
functions. 


2.  Boundary  Collocation 

The  boundary  collocation  technique  developed  by  Gluckman,  Pfeffer  and  Weinbaum 
(1971)  has  been  used  to  solve  a  wide  variety  of  low-Reynolds-number  problems 
where  the  system  boundaries  do  not  conform  to  a  single  orthogonal  co-ordinate 
system.  A  related  technique  was  developed  by  O'Brien  (1968)  to  calculate  the 
flow  past  a  slightly  deformed  sphere.  The  earliest  applications  of  the 
technique  were  limited  to  axisymraetric  problems  which  were  solved  using  the 
stream  function.  Since  then,  the  technique  has  been  applied  directly  to  the 
Stokes  equation  and  three  dimensional  problems  including  the  sedimentation  of 
three  spheres  with  centers  in  a  vertical  plane  (Ganatos,  Pfeffer  and  Weinbaum 
1978),  the  motion  of  a  sphere  between  two  parallel  Infinite  plates  (Ganatos, 
Pfeffer  and  Weinbaum  1980)  and  the  sedimentation  of  a  sphere  in  an  inclined 
channel  (Ganatos,  Weinbaum  and  Pfeffer  1982).  For  our  two-sphere  problem,  a 
suitable  co-ordinate  system  exists.  Nevertheless  we  use  the  collocation 
technique  because  accurate  numerical  results  are  obtained  with  minimal  human 
computation. 

The  essential  idea  behind  the  collocation  technique  is  as  follows.  The 
velocity  field  can  be  represented  by  an  expansion  in  terms  of  basis  functions, 
each  of  which  satisfies  the  equations  of  motion.  In  general,  the  number  of 
elements  in  the  basis  set  is  not  finite  because  of  the  interactions  between 
the  spheres.  However,  the  higher  order  elements  are  usually  unimportant. 
Consequently,  the  series  can  be  truncated  at  N  terras  and  the  coefficients  of 
the  retained  basis  functions  determined  by  setting  the  boundary  condition  at  N 
collocation  points  (hence  the  name  boundary  collocation).  It  is  found 
empirically  that  the  lower  order  coefficients  converge  rapidly  as  N  is 
increased.  This  is  important  since,  as  shown  later,  the  force,  torque  and 
stresslet  depend  only  on  the  first  and  second  order  coefficients. 
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2.1  The  Velocity  Representation 

The  velocity  field,  y,  satisfies  the  Stokes  equation: 

-Vp  +  yV2 v  ■  0 

and  the  equation  of  continuity: 

V*v  -  0. 


(2.1) 


(2.2) 


The  boundary  conditions  are  those  associated  with  each  resistance  problem. 

Our  goal  in  this  section  is  to  construct  general  forms  of  the  velocity 
representations  and  boundary  conditions  which  together  encompass  the  complete 
set  of  resistance  problems.  Then,  any  resistance  problem  of  interest  can  be 
obtained  by  selecting  the  appropriate  set  of  parameters.  This  structure  is 
readily  passed  on  to  the  computer  codes  and  the  result  is  a  versatile,  yet 
short  subroutine  (less  than  60  lines  of  code)  which  calculates  the  entire 
collection  of  functions. 

As  shown  in  Happel  and  Brenner  (1965),  the  disturbance  velocity  field  can 
be  represented  using  Lamb's  general  solution: 


v(x)  -  v°°(x)  -  I  {  V*_n_1  +  Vx(xx_n_1 ) 
n»1 


(2.3) 


(n+1 ) 
n(2n-1 ) 


*  P-n-1  ‘ 


n-1  2n(2n-1) 


—  I  • 


where  P.n_1 t  *_n_i  and  X_n_1  ar®  exterior  spherical  harmonics. 

Following  Ganatos  et.  al.  (1978),  the  velocity  is  written  as  a 
superposition  of  two  expansions,  one  centered  at  sphere  1  and  the  other  at 
sphere  2. 
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The  spherical  harmonics  are  expanded  as: 
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For  each  resistance  function,  we  will  actually  require  only  one 
particular  m,  the  value  which  appears  in  the  surface  velocity,  and  the 
^-dependence  will  factor  from  our  problem  as  shown  in  the  next  subsection. 
Thus  for  two-sphere  problems,  a  one-dimensional  collocation  (in  6)  is 
possible,  even  if  the  flow  is  three-dimensional. 


2.2  Application  of  the  Boundary  Conditions 

The  disturbance  field  must  decay  far  away  from  both  particles.  At  the  sphere 
surface,  the  disturbance  velocity  must  equal  a  surface  velocity,  v  ,  which  is 
the  difference  between  the  particle's  rigid-body  motion  and  the  ambient 
velocity.  Thus  all  relevant  cases  are  included  in  the  following  expression  of 
the  boundary  condition  at  sphere  1: 

2  £ 

*s  '  1-1  m'o  [  ',|rlP?<oos9)CiWo.  *  Atn,slnm*]  1  (2‘5) 

+  7*{C-|  r^P®(cose)Bimcosm$]  }  ]. 


£*7 


,  «*4  .*  »V  **„  .,*«  -  .  -  .  *«  *  .%/.  *T.  *  »  ,  r.  ^  p  »  »\  *,  **,  ».  *,  «  _ 

W>  J *k Xj, Vj|  »,»  _>  -N  j,  -  .A  _NjA  /*,*•  .  -  A-Vu'.-N  /•  ,  *  ,>  ,**  ,%  ,  -  /-O 


In  (2.5)  the  cosm<j>  terms  are  omitted  for  mSI ,  since  those  resistance  problems 
are  equivalent  to  those  obtained  from  the  sinm$  terms  (with  a  rotational 
co-ordinate  transformation  about  the  sphere-sphere  axis).  The  following  table 
shows  the  required  v  for  each  resistance  function. 
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Table  1 . 


Non-zero 

Coefficient 


Resistance 

Function(s) 


Translation  along  sphere-sphere  axis. 
Translation  perpendicular  to  axis. 
Axisymmetric  straining. 


Rate-of-strain  as  in  ZX  shear  flow. 


Hyperbolic  straining  in  XY  plane. 

Rotation  about  sphere-sphere  axis.  . 

Rotation  with  axis  perpendicular  to 
sphere-sphere  axis. 


YaB  YaB 
XA  ,XG 

vaB  va8  yoB 
*A  •  XB  ’  *G 

xJX2 


vaB  vaB 
XC  ■  XH 


We  now  examine  the  form  taken  by  Lamb's  representation  at  the  sphere 
boundaries.  We  use  the  cylindrical  coordinate  system  (z,R,+)  as  shown  in 
Figure  1.  The  z,  R  and  $  velocity  components  in  equation  (2.4)  are  equated  to 
the  corresponding  components  of  the  surface  velocity  in  (2.5).  Dependence  on 
<f>  occurs  for  problems  with  m21 ,  but  factors  out  as  follows:  a  factor  of  sinm$ 
in  the  z-component  and  R-component  equations  and  a  factor  of  cosmt  in  the 
^-component  equation.  Thus  the  boundary  conditions  on  sphere  1  (i.e.  r1  -  1) 
requires  that: 
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(2.6c) 


where  5a  ■  C03eft.  There  is  an  analogous  set  of  equations  from  sphere  2. 
Equations  (2.6a)  and  (2.6c)  follow  directly  from  the  z-component  and 
^-component  equations.  Equation  (2.6b)  is  obtained  by  subtracting  the 
^-component  equation  from  the  R-coraponent  equation  This  last  maneuver  allows 
us  to  remove  the  singularities  from  the  poles.  Further  details  are  given  in 


Section  2.5. 
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The  series  is  now  truncated  at  N  terms.  The  6N  unknown  coefficients, 
112  2  2 

b  ,  c  ,  a  ,  b  and  c  ,  are  determined  by  applying  the  truncated 
mn  mn  mn  mn  mn 


version  of  (2.6a-c)  at  2N  collocation  points  on  the  surfaces  and  solving  the 
resulting  6Nx6N  system.  We  remind  the  reader  that  the  parameters  m  and  £  are 
specified  by  the  resistance  problem.  In  the  next  section,  the  computations 
are  simplified  by  exploiting  the  mirror  symmetry  with  respect  to  the  XY  plane. 

2.3  Mirror  Symmetry  about  the  XY  Plane 

For  the  general  problem  of  two  unequal  spheres,  it  has  been  shown  that  the 
larger  sphere  requires  more  points  (Liao  and  Krueger  1980).  For  two  equal 
spheres,  the  points  are  distributed  in  equal  numbers  between  the  two,  at 
equidistant  spacings  (Gluckman,  Weinbaum  and  Pfeffer  1971).  Furthermore,  we 
decompose  each  resistance  problem  into  subproblems  that  exploit  the  fore-aft 
mirror  symmetry  with  respect  to  the  XY  plane.  Symmetry  dictates  that  the 
coefficients  in  the  series  centered  at  sphere  1  either  equal  or  are  negatives 
of  the  corresponding  coefficients  in  the  other  series.  This  also  holds  for 
the  truncated  expansion  as  long  as  the  collocation  points  on  sphere  2  are 
placed  at  the  mirror  images  of  the  points  chosen  for  sphere  1 . 

An  examination  of  the  resistance  problems  reveals  that  they  either  posess 
one  of  the  following  two  types  of  symmetry  or  may  be  decomposed  into  two 
subproblems,  with  a  subproblem  of  each  symmetry  type.  A  velocity  field  with 
mirror  symmetry  with  respect  to  the  XY  plane  satisfies: 


vx(x,y,z)  -  vx(x,y,-z), 
Vy(x, y, z)  -  vy(x,y,-z), 
vz(x,y,z)  -  -vz(x,y,-z). 


(The  flow  vectors  in  the  half-spaces  are  mirror  images  of  each  other).  A 
field  with  mirror  anti-symmetry  satisfies: 


V 


vx(x,y,z)  -  -vx(x,y,-z), 
vy(x,y,z)  -  -vy(x,y,-z), 
v2(x,y,z)  -  v2(x,y,-z). 


For  problems  with  these  symmetries,  the  coefficients  for.  the  terms 
centered  at  sphere  2  are  given  by: 


where  the  symmetry  parameter,  S,  is  defined  by: 


S  -  { 


1  for  problems  with  mirror  symmetry, 

-1  for  problems  with  mirror  anti-symmetry. 


and  the  collocation  equation  from  a  point  on  the  surface  of  sphere  2  becomes 
identical  to  that  from  the  image  point  on  sphere  1.  Thus  equations  (2.6a-c) 
and  their  counterparts  from  sphere  2  reduce  to: 
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Thus  each  resistance  function  is  calculated  by  solving  an  appropriate  set  of 
problems  of  the  form  given  by  (2.7a-c).  We  shall  see  that  each  resistance 
function  is  a  linear  combination  of  a  small  number  of  the  lower  order 
coefficients. 

The  remainder  of  Section  2  contains  information  relating  to  the  code 
development  from  (2.7a-c).  Readers  who  are  more  interested  in  the  end 
applications  may  prefer  to  skip  to  Section  3. 

2.4  Axlsymraetrlc  Problems 

In  axisymmetric  problems  certain  velocity  components  vanish  identically 

because  of  symmetry.  The  3Nx3N  collocation  system  reduces  to  either  an  N*N  or 

2 N*2N  system.  Axisymmetric  problems  with  v  -0  have  m-0  so  that  equation 
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(2.7c)  and  c  vanish  identically  leaving  a  2Nx2N  system.  On  the  other  hand, 
mn 

axisymmetric  swirl  problems  with  v^  as  the  only  nonvanishing  component  reduce 

v 

to  an  N*N  system  since  a  ,b  ,  (2.7a)  and  (2.7b)  vanish  identically.  The 

mn  mn 

reduced  system  can  be  obtained  from  the  c^  terms  of  (2.7b)  (i.e.  the  2-3 

mn 

block  of  the  larger  system). 

2.5  Stability  of  the  Collocation  System  of  Equations 

In  this  section,  we  review  earlier  reports  on  the  stability  of  the  collocation 
system  of  equations  and  present  new  findings. 

Problems  occur  in  the  collocation  scheme  if  points  are  placed  at  the 
equator.  el  -  ir/2,  and  the  poles.  If  a  point  is  placed  at  the  equator,  the 
terms  from  sphere  1 ,  which  are  normally  greater  than  those  from  sphere  2 
vanish  (since  cose^  vanishes  at  the  equator)  thus  destabilizing  the  system. 
This  problem  was  circumvented  in  the  original  work  by  Gluckman  et.  al.  (1971) 
by  using  twin  points  at  89°  and  91°,  an  approach  which  they  justified  by 
examining  the  convergence  behavior,  in  the  limit  of  small  e,  for  twin  points 
at  90-e  and  90+e  degrees. 

At  the  poles,  the  system  is  indeterminate  because  equations  from  the  R 
and  $  components  become  identical.  Gluckman  et.  al.  (1971)  avoided  this 
problem  by  not  placing  any  points  at  the  poles. 

We  have  avoided  the  problem  at  the  equator  by  using  an  even  number  of 
points  spaced  at  equidistant  intervals.  The  Indeterminacy  at  the  poles  was 
removed  by  taking  the  difference  of  the  R-component  and  ^-component  equations 
to  arrive  at  (2.6b)  and  (2.7b).  The  source  of  the  indeterminacy  is  then 
apparent.  Equation  (2.7a-c)  have  zeros  of  multiplicity  m,  ra+1  and  m-1  at  the 
poles.  They  may  be  removed  by  factoring  sin“e1 ,  sinn+1 e1  and  sinm-1e 


respectively. 


When  the  spheres  are  far  apart,  accurate  solutions  are  obtained  with  as 
few  as  four  points.  Tests  show  that  our  collocation  scheme  compares  favorably 
with  that  of  Gluckman  et.  al.  (1971)  for  R/a  between  2.1  and  10.0. 

Furthermore,  in  the  strong  interaction  region,  our  scheme  converges  faster 
because  the  error  profile  in  the  gap  region  is  reduced  by  the  "Hermite 
interpolation"  nature  of  the  approxiraant.  Although  we  do  not  place  a  point  at 
the  equator,  the  large  number  of  points  used  in  the  near-field,  e.g.  N-60, 
ensures  the  presence  of  collocation  points  near  the  equator. 

The  code  development  required  a  system  solver  and  a  routine  for  the 
spherical  harmonics.  We  used  the  LINPAK  routines  DGECO  and  DGESL  to  invert 
the  system.  However,  essentially  identical  results  (15  significant  figures) 
were  obtained  with  other  (slower)  routines.  A  recursion  scheme  was  developed 
for  the  spherical  harmonics  because  our  application  required  the  harmonics 
divided  by  factors  of  sinme.  The  stability  of  these  routines  was  spot  checked 
by  comparison  with  the  tables  in  Abramowitz  and  Stegun  (1964). 


3-  Calculation  of  Resistance  Functions 

In  this  section,  we  extract  the  resistance  functions  from  the  Information 
contained  in  the  collocation  solution.  There  are  several  methods  for 
obtaining  the  force,  torque  and  stresslet  on  a  particle  from  the  velocity 
solution  (see  Chapter  3  of  Happel  and  Brenner,  1965).  We  present  here  a 
simple  but  powerful  method. 

The  velocity  field  which  was  previously  represented  by  Lamb’s  general 
solution  can  also  be  represented  by  the  twin  multipole  expansion  (Jeffrey, 
1974): 

2 

v  -  l  1  F  -  (S  +  T  )  *V  ♦  •••  }  »I(x-x  )/(8iru) , 
ot  *oi  -a  •  o 

a-1 

with  the  Oseen  tensor,  I  defined  by 
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This  representation  is  useful  because  it  reveals  that  the  force  on  particle  a 

appears  as  the  coefficient  in  the  term  that  decays  as  |x_?a|  1  while  the 

dipole  moments,  and  -  -  ,^S*Ia  appear  as  the  coefficients  in  the  terms 

•2 

that  decay  as  | x— x^j  .  More  explicitly,  in  the  notation  of  Chwang  and  Wu 
(1975),  the  force,  torque  and  stresslet  appear  as  the  coefficient  of  the 
Stokeslet,  Rotlet  and  Stresslet.  Therefore,  we  rearrange  the  terms  in  Lamb's 
representation  to  form  the  Stokeslet,  Rotlet  and  Stresslet  and  obtain  the 
relation  between  the  force  and  dipole  moments  and  the  coefficients  (in  Lamb's 
representation). 

The  above  ideas  are  put  into  practice  for  each  resistance  problem  in  the 
following  two-step  procedure. 

1)  In  the  first  step,  we  take  the  arbitrary  but  convenient  convention  of 
setting  the  appropriate  A,  or  B,  equal  to  one  (and  set  all  others 
equal  to  zero). 


V 


2)  The  force,  torque  and  stresslet  in  the  multipole  expansion  are  expressed  in 
terms  of  the  resistance  functions  as  given  by  (1.5a-f).  This  expansion  has 
the  same  form  as  the  (rearranged)  Lamb's  representation. 

After  the  above  prescribed  algebra,  we  arrive  at  formulae  for  the  resistance 

functions  in  terms  of  the  coefficients  in  Lamb's  solution  -  a  (i,,m,S), 

mn 

b  (l,m,S)  and  c  (l,m,S).  The  arguments,  £,,  m  and  S  indicate  which 
mn  mn 

coefficient  is  retained  in  the  surface  velocity  (see  Table  1  on  page  13)  and 
the  type  of  mirror  symmetry.  As  discussed  earlier,  each  function  requires  a 
superposition  of  a  mirror  symmetric  and  mirror  anti-symraetric  problem,  with 
the  exception  of  the  scalar  functions  from  M  which  already  podess  the 
symmetry.  As  shown  below  the  mirror-symmetric  solutions  are  summed  to  give 
the  1-1  functions  and  their  difference  is  taken  to  give  the  1-2  functions. 

The  functions  associated  with  translational  motions  and  rate-of-strain 
fields,  with  the  argument  (S.,m,S)  denoting  which  A^^  is  set  equal  to  one,  are 


given  by: 
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(3.3a) 


(3-3b) 


(3.30 


(3.3d) 


(3-3e) 


(33f ) 


(3.3g) 


(3- 3h) 


(3.3D 
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Y°2  -  ^Ca12(1t1.-1)  -  a12(1,1,1)] 


*?t  *  ^2  ■T^02<2-°-1) 


■  -  a  (2  1  -1 ) 

X11  *12  To®^2’1’  V 


7M  ♦  -  - q  (22  1 ) 

Z11  Z12  T0®22U,2,.J* 


(3.3J) 

(3.3k) 

(3-30 

(3.3m) 


The  functions  associated  with  rotational  motions,  with  the  argument  (£,m,S) 
denoting  which  is  set  equal  to  one,  are  given  by: 


*n  "  5Cc01(1,0,"1)  +  coi(1'°*1)]  (3.3n) 

xf2  “  -  ^tc01(1,0,-1)  -  c01 (1 .0,1 ) 3  (3.3o) 

Y^  -  iCcn(1,1,-1)  ♦c11(1,1,1)]  (3.3p) 

*12  "  2C°11<1,1»“1>  ~  e11(1»1*1)3  (3.3q) 

Y^i  -  -  Jta12(1(i;-1)  +  a12(1,1,1)]  (33r) 

Y^2  “  *  ^Ca12(1,1,-1)  -  a12(1,1,1)].  (3.3s) 


This  completes  the  calculation  of  the  resistance  functions.  The  mobility 
functions  are  calculated  from  the  resistance  functions  as  prescribed  by 
equations  (1.10)  to  1.15). 
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4.  Results  and  Conclusions 


In  this  section,  we  examine  the  results  for  the  resistance  and  mobility 


functions  and  compare  them,  where  available,  with  those  obtained  by  other 


means.  The  resistance  and  mobility  functions  are  plotted  vs.  the  reduced 


sphere-sphere  separation,  R/a,  in  Figures  2a-fc  and  3a-£,  respectively.  These 


plots  were  generated  using  12  collocation  points  for  R£3  and  24  points  for 


R<3.  At  these  levels,  convergence  has  been  obtained  far  beyond  the  resolution 


of  the  plotter.  We  note  that  ,  Y^+Y^2,  y°2  and  y°  are  not  monotonic 


functions  but  have  extrema.  This  phenomenon  is  consistent  with  the 


requirement  that  two  neutrally-buoyant  spheres  move  as  a  rigid  body  when  in 


contact. 


For  comparative  purposes,  an  extensive  table  (not  presented  here  but 


available  upon  request  from  SK)  was  created  for  all  functions.  In  the 


construction  of  the  table,  the  number  of  collocation  points  was  increased 


until  convergence  was  obtained  to  five  significant  figures.  This  table  was 


used  as  "the  collocation  results"  in  the  comparison.  Except  for  almost 


touching  spheres,  the  number  of  collocation  points  and  therefore  the  order  of 


the  system  of  equations  was  well  within  the  memory  limitations  of  a  VAX 


11/780.  At  sphere-sphere  gaps  of  0.01a  the  most  difficult  function,  , 


converged  to  three  significant  figures  at  60  collocation  points.  The  program 


in  its  present  form  allows  up  to  100  points. 


In  the  far-field,  asymptotic  solutions  are  either  available  or  readily 


obtained  by  the  method  of  reflections.  The  collocation  results  matched  these 


far-field  solutions  for  all  functions.  In  fact,  the  far-field  solutions  were 


the  primary  defense  against  program  bugs. 


A  more  stringent  test  was  available  for  resistance  and  mobility  functions 


from  A,  B,  C,  a,  b  and  c  because  of  the  recent  twin-multi  pole  expansion 


xl la  AND  yl  la  3b.  x!2a  AND  y!2a 
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results  of  Jeffrey  and  Onlshl  (1984).  The  collocation  results  matched  those 
obtained  from  a  fifteen-term  form  of  their  expansion  solution.  The  agreement 
was  exact  (in  the  sense  that  it  was  limited  only  by  machine  roundoff  errors 
e.g.  over  12  significant  figures)  except  at  the  small  separations  mentioned 
above,  where  both  techniques  require  special  modifications. 

Accurate  results  are  also  available  for  certain  combinations  of  the  g 
functions.  The  relative  velocity  between  two  spheres  in  a  shear  field  which 
was  obtained  by  Lin,  Lee  and  Sather  (1970)  using  bispherical  coordinates 
furnishes  a  test  for  x®j-x®2  and  y^-j^y^*  Here  a8ain»  agreement  was  obtained 
to  all  significant  figures  presented  in  the  earlier  work. 

In  the  near-field,  the  collocation  results  for  A,  B,  C,  a,  b,  and  c  were 
compared  with  the  near-field  lubrication  solutions.  The  latter  are  collected 
from  the  literature  and  presented  in  Jeffrey  and  Onishi  (1984).  Plots  for  the 
mobility  functions  presented  in  Figures  4a-j  show  that  the  two  solutions  match 
in  an  overlap  region,  but  also  show  that  the  collocation  technique,  at  least 
in  its  current  form,  cannot  "turn  the  sharp  corner"  in  the  y  functions. 

In  conclusion,  we  believe  that  we  have  successfully  calculated,  using  the 
boundary  collocation  technique,  the  complete  set  of  resistance  and  mobility 
functions  required  for  determining  the  force,  torque,  stresslet  on  and  motions 
of  two  equal  sized  spheres  in  an  ambient  velocity  composed  of  a  uniform  stream 
and  linear  field.  The  results  are  accurate  over  all  separations  except  at 
almost-touching  (R/a  <2.01).  In  particular,  this  report  provides  an  accurate 
algorithm  for  the  computation  of  the  stresslet  functions  over  a  wide  range  of 
sphere-sphere  separations. 
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5.  Sample  Applications 


5.1  The  Bulk  Stress  In  a  Suspension  of  Spheres  to  Order  o 

Batchelor  and  Green  (1972b)  have  derived  the  following  rigorous  expression  for 
the  viscosity  of  a  suspension  of  Identical  rigid  spheres  in  a  steady  pure 


straining  motion. 


W  -  •  *  h  *  =2l!  * 


■f s\  J(c)q(c)c2<ic  }  . 


(5.1) 


where  c  »  R/a  and  J(()  and  q(;)  are  determined  from  the  mobility  functions  as 


J(C)  -  >1  -  x“  ♦  2y”  ♦  2z?  }  , 


(5.2) 


q(t)  -  [1  -  A(c)]-1  exp{/c  jjffi-rjfldx}. 


(5.3) 


A(c)  and  B(c)  are  mobility  functions  that  arise  in  the  expression  for  the 
relative  velocity  between  the  two  spheres  and  are  as  follows: 


A(?)  -  4(x®1-x®2)/;, 


(5.4) 


BU)  -  8(y8ryf2)/c 


(5.5) 


At  the  time  of  their  work,  information  on  J  was  limited  to  the  far-field 
region  (obtainable  by  the  method  of  reflections)  and  the  value  at  touching. 

As  mentioned  in  their  paper,  the  interpolation  of  the  J  curve  in  the  region 
2. 0025<c<3  was  the  primary  source  of  uncertainty  in  the  final  numerical 
result. 

Figures  5a-5d  show  plots  of  J  and  also  their  stresslet  functions  K,L  and 
M.  The  dashed  line  is  the  result  obtained  using  the  method  of  reflections. 
The  nxn  at  R/a  -  2  indicates  the  values  for  touching  spheres.  As  is  the  case 
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Figure  5.  Comparison  of  the  collocation  and  method  of  reflection  solutions 
for  the  stresslet  functions  of  Batchelor  and  Green. 
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with  many  of  the  mobility  problems,  we  see  that  the  method  of  reflections 
result  is  accurate  to  quite  small  separations.  The  explanation  is  that  in 
mobility  problems,  the  leading  terms  in  the  higher  order  reflections  are  from 
the  relatively  weak  dipole-dipole  interactions.  We  note  that  L  and  M  are  not 
monotonlc. 

When  (5.1)  is  evaluated  with  the  collocation  result,  we  find  that  the 
integral  is  less  than  the  earlier  estimate.  Specifically,  the  integral  from 
2.0025  to  3.0  is  found  to  be  0.384  instead  of  0.449.  Thus  the  coefficient  of 
the  c  term  is  7-1  instead  of  7.6,  i.e., 

Ueff/u  -  1  *  fc  +  7.1c2. 

The  new  result  is  within  the  ±10$  uncertainty  bound  stated  by  Batchelor  and 
Green. 

5.2  The  motion  of  a_  rigid  dumbbell  in  a  shear-field 

Rigid  dumbbells  have  been  used  to  model  suspensions  of  polymers  with  stiff 
backbones  (see  Bird  et.  al.  1977).  In  their  models  the  hydrodynamic  problem 
is  simplified  by  assuming  that  the  dumbbell  consists  of  two  point  forces 
connected  by  a  rigid  rod.  (The  rod  has  no  hydrodynamic  resistance).  In  this 
section,  we  show  that  if  one  replaces  the  point  forces  with  spheres,  the 
necessary  hydrodynamic  functions  can  be  extracted  from  our  two-sphere 
functions. 

The  motion  of  a  neutrally  buoyant,  axisymmetric  particle  in  a  shear  field 
was  completely  solved  by  Bretherton  (1962)  who  showed  that  the  motion  of 
almost  all  axisymmetric  particles  was  the  same  as  that  of  some  ellipsoid  of 
revolution.  Thus,  for  our  dumbbells,  we  use  our  two-sphere  functions  to  find 
the  "effective  spheroid". 


The  geometry  of  a  symmetric  dumbbell  is  completely  determined  by  the 
center-to-center  separation,  R  and  sphere  radii,  a.  Furthermore,  if  all 
distances  in  the  problem  are  scaled  with  "a",  then  the  possible  dumbbell 
shapes  are  spanned  by  varying  R/a  over  [2,»).  If  the  dumbbell  rotates  about 
its  center  of  mass  as  «xx,  then  the  spheres  move  as: 

U  -  U  +  w*x  ,  a-1,2,  (5.6) 

-a  -  -  -a 

We  can  now  write  the  forces  and  torques  on  each  sphere  in  terms  of  the 
two-sphere  resistance  functions.  The  torque  on  the  dumbbell  is  the  sum  of  the 
torques  on  each  sphere  plus  the  couple  from  the  forces  on  the  spheres.  After 
some  algebra,  we  arrive  at: 

Ti  "  Cij(“j‘V  ‘  EijkEJk’  (5'7) 

where  the  dumbbell  resistance  functions  can  be  written  in  terms  of  the 
two-sphere  resistance  functions  as: 

°ij  •  *CVj  *  lC(su  -  Vj>-  (5-8) 

with  XC  -  2(X^+xf2) 

and  YC  -  -1r2(Y^-Y*2)  +  2R(Y®  -Y®2)  +  2<yJ1+Y^2). 

Hljk  "  Y  (djeki£+dkeJU)d£.  (5.9) 


YH  -  ^(Y^-Y^)  ♦  ^(Y^-Y^)  -  R(Y°-Y°2)  +  2(Y“1+Y“2), 


wi  th 


If  we  set  T  -  0  in  (5.7),  we  can  solve  for  the  angular  velocity  as: 

a;  -  Q  ♦  {2YH/YC}dK(E*d).  (5.10) 

A  prolate  spheroid  with  aspect  ratio  p  rotates  as  given  by  (5.10),  but  with 
2YH/YC  replaced  by  (p2-1 )/(p2+1 ) .  Therefore,  the  dumbbell  moves  like  a 
prolate  spheroid  with  the  aspect  ratio, 

p  -  [(1  +  2YH/YC)/(1  -  2YH/YC)]1/2.  (5.11) 

Plots  of  the  "equivalent-spheroid"  aspect  ratio  vs.  the  dumbbell  shape 
parameter,  R/a  are  shown  in  Figure  6. 
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Appendix 


1.  Notes  on  the  Computer  Programs 

The  following  section  contains  listings  of  five  programs.  The  "main"  program, 
MANDR. FOR  requires  only  two  inputs:  the  number  of  collocation  points  and  the 
scaled  separation  between  the  spheres,  R/a.  MANDR. FOR  occupies  five  pages  and 
this  length  may  give  the  wrong  impression  that  the  algorithm  is  Involved. 
However  this  program  merely  sets  up  the  appropriate  inputs  for  the  entire 
collection  of  resistance  functions  and  passes  them  to  subroutine  STOKES. FOR. 
The  heart  of  the  algorithm  is  contained  in  STOKES. FOR,  which  occupies  only  a 
single  page.  If  one  were  interested  in  a  particular  resistance  function,  the 
correct  set  of  input  parameters  for  STOKES. FOR  may  be  deduced  by  examining  the 
appropriate  module  in  the  main  program. 

As  mentioned  in  Section  2,  axisymmetrio  problems  reduce  to  smaller 
systems.  In  axisymmetric  translational  and  straining  problems,  the  resulting 
2Nx2N  system  consists  of  the  upper-left  blocks  of  the  general  system. 
Consequently,  in  STOKES. FOR  we  need  only  insert  the  line: 

IF  (M  .EQ.  0)  IDIM  -  2*NPT 

before  calling  the  inversion  (LINPAK)  subroutines. 

In  axisymmetric  swirl  problems,  the  N*N  system  is  not  located  in  the 
upper-left  portion  of  the  general  system.  Consequently,  the  appropriate  block 
cannot  be  obtained  as  easily  from  the  general  case.  We  have  taken  the  simple 
remedy  of  taking  the  2-3  block  out  of  STOKES. FOR  and  creating  the  Regenerate 
subroutine,  STOKESD.FOR. 

The  following  parameters  are  passed  in  and  out  of  STOKES. FOR: 

R:  Sphere-sphere  separation  (scaled  by  sphere  radius). 


/  J 


NSTAR:  Order  of  the  lowest  non-zero  multipole.  (e.g.  if  the  hydrodynamic 
force  on  the  sphere  is  nonzero,  NSTAR-1 .  If  the  force  is  zero, 
then  NSTAR  is  2  in  our  application  because  of  the  dipole  moment. 

M:  Parameter  m  from  the  resistance  problem. 

SGN:  Symmetry  parameter  S. 

NPT:  Number  of  collocation  points  (N  in  the  text). 

RHS:  One  dimensional  array  which  contains  the  RKS  of  the  system  of 

equations.  It  originates  from  the  surface  velocity  evaluated  at 
the  N  collocation  points. 

AMN:  One  dimensional  output  array  of  the  a  (with  n  *  NSTAR,  ...  ). 

mn 

BMN:  One  dimensional  output  array  of  the  b  (with  n  -  NSTAR,  ...  ). 

mn 

CMN:  One  dimensional  output  array  of  the  c  (with  n  •  NSTAR,  ...  ). 

mn 

Function  subroutine  PNS(N,M,X)  calculates  Pm(X)[sine]~m.  (PNS  stands  for 

n 

P-No-Sine).  PNSP(N,M,X)  does  the  corresponding  steps  for  p“'(X).  Subroutine 
DMATIN.FOR  is  used  when  the  mobility  functions  are  calculated  from  the  inverse 
of  the  resistance  matrix. 

If  these  programs  are  installed  on  another  computer,  the  sample  run 
provided  after  the  program  listings  may  be  used  to  test  the  installation. 
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MOBILITY  AND  RESISTANCE  FUNCTIONS  FOR  2  IDENTICAL  SPHERES 
BY  BOUNDARY  COLLOCATION  OF  LAMB'S  GENERAL  SOLUTION 

SANGTAE  KIM 

DEPARTMENT  OF  CHEMICAL  ENGINEERING 
AND  MATHEMATICS  RESEARCH  CENTER 
UNIVERSITY  OF  WISCONSIN 

RICHARD  T.  MIFFLIN 
DEPARTMENT  OF  CHEMICAL  ENGINEERING 
PRINCETON  UNIVERSITY 

VERSION  2:  APRIL  24.  1984 


PROGRAM  WAS  DEVELOPED  AND  TESTED  ON  A  VAX  11/780.  1985-1984 

MAIN  PROGRAM:  MANOR 

SUBROUTINES  CALLED  FROM  MAIN:  STOKES,  STOKESD,  DMATIN 
SUBROUTINES  CALLED  FROM  STOKES.  STOKESD:  LINPAK  ROUTINES  DGECO.DGESL 

PNS.PNSP 


PROGRAM  DESCRIPTION 


MAIN  PROGRAM  SETS  PARAMETERS  CORRESPONDING  TO  AMBIENT  VELOCITY  FIELD 
FOR  EACH  RESISTANCE  FUNCTION.  AFTER  CALCULATION  OF  THE  MULTIPOLES. 
MAIN  ALSO  SCALES  THE  RESULT  ACCORDING  TO  THE  NON-OIMENSIONALIZATION 
FOR  EACH  RESISTANCE  FUNCTION.  MAIN  CALCULATES  THE  MOBILITY  FUNCTION 
BY  INVERTING  THE  GRAND  RESISTANCE  MATRIX. 

GIVEN  THE  PARAMETERS  FOR  THE  AMBIENT  VELOCITY,  SUBROUTINE  STOKES 
RETURNS  MULTIPOLES  COEFFICIENTS.  THE  SYSTEM  OF  ECU AT IONS  IN  THIS 
SUBROUTINE  ARE  OBTAINED  BY  APPLYING  THE  BOUNDARY  CONDITIONS  AT  EACH 
COLLOCATION  POINT. 

SUBROUTINE  STOKESD  IS  A  SUBSET  OF  STOKES  AND  IS  USED  FOR  THE 
DEGENERATE  CASES  INVOLVING  THE  RESISTANCE  FUNCTIONS  XI 1C  AND  X12C. 

PNS  AND  PNSP  ARE  DERIVED  FROM  THE  ASSOCIATED  LEGENDRE  FUNCTIONS 

IMPLICIT  DOUBLE  PRECIS I ON(A-H.O-Z) 

DIMENSION  RHS(300) ,AMN(100) ,BMN(100) ,CMN(100) , 

0  COS1 (100) ,AY(4,4) , DUM4Y(4, 1 ) 


READ  IN: 

1)  NUMBER  OF  COLLOCATION  POINTS  (EVEN  INTEGER):  IPTS 

2)  CENTER  TO  CENTER  SEPARATION  BETWEEN  SPHERES:  R 

READ  (40.5)  IPTS.R 
5  FORMAT  (I10.F10.0) 

SET  COLLOCATION  POINTS: 

PI  -  3. 141S9265358979D0 
DIPTS  -  DFLOAT(IPTS-I) 

DO  10  1-1. IPTS 
THETA  -  DFLOAT ( 1—1 )/DIPTS»PI 
10  COSI(I)  -  DCOS(THETA) 


o  o  o  o  o  o  o  o 


•  MM....  CALCULATION  OF  XI 1A.  X12A.  X11C.  X12C  •• 

DO  20  I-1.IPTS 

12  -  I+IPTS 

13  -  I2+IPTS 
RHS(I)  -  1 .DO 
RHS(I2)  -  0.D0 

20  RHS(I3)  -  0.D0 

CALL  ST0KES(R, 1 .0,-1 . DO, IPTS.RHS, AMN.BMN.CVM) 
T1  -  2.DO/3.DO*AMN(1) 

T3  -  -.5D0*AMN(2) 

DO  30  1-1.IPTS 

12  -  I+IPTS 

13  -  I2+IPTS 
RHS(I)  -  1 .DO 
RHS(I2)  -  0.D0 

30  RHS(I3)  -  0.D0 

CALL  ST0KES(R, 1.0,1. DO, IPTS.RHS.AMN.BMN.CMN) 
T2  -  2.D0/3.D0*AMN(1 ) 

T4  -  -.5DO»AMN(2) 

XI 1A  -  .SD0*(T1  +  T2) 

X12A  -  .5D0*(T1  -  T2) 

XI 1G  -  .500* (T3  +  T4) 

XI 2G  -  .500* (T3  -  T4) 


.  CALCULATION  OF  Y11A,  Y12A.  Y11B.  Y12B.  Y11C,  Y12G 

DO  40  I-1.IPTS 

12  -  I+IPTS 

13  -  I2+IPTS 
RHS(I)  -  0 .  DO 
RHS(I2)  -  0.D0 

40  RHS(I3)  -  1.00 

CALL  ST0KES(R.1, 1.-1. DO. IPTS.RHS.AMN.BMN.CMN) 

T1  -  2.DO/3.DO*AMN(l) 

T3  -  2.D0*CMN(1) 

T5  -  -.5D0*AMN(2) 

DO  SO  I-1.IPTS 

12  -  I+IPTS 

13  -  I2+IPTS 
RHS(I)  -  O.DO 
RHS (12)  -  O.DO 

50  RHS(I3)  -  1 .DO 

CALL  STOKES(R. 1.1,1. DO. IPTS. RHS. AMN.BMN.OM) 

T2  -  2.D0/3.D0*AMN(1) 

T4  -  2.00*CMN(1) 

T6  -  -.5D0*AMN(2) 

Y11A  -  .5D0*(T1  +  T2)  , 

Y12A  -  .500. (T2  -  T1 ) 

Y11B  -  -.5D0»(T3  +  T4) 

Y12B  -  . 5D0»(T3  -  T4) 

Y11C  -  .5D0»(T5  +  TO) 

Y12C  -  . 5D0»(T6  -  T5) 


no  o  o  o  o  o  o  oooo  ooo 


C 


00  120  I-1.IPTS 

12  -  I+IPTS 

13  -  I2+IPTS 
RHS(I)  -  0.00 
RHS (12)  -  0.00 

120  RHS(I3)  -  6.00 

CALL  ST0KES(R,2.2.1.D0,IPTS,RHS.AMN,BMN.CMN) 

Z1112M  -  . 1D0*AMN(1 ) 

■  '  END  OF  CALCULATIONS  FOR  THE  RESISTANCE  FUNCTIONS  - 

WRITE  (41,500)  IPTS.R 
WRITE  (41.510) 

WRITE  (41,520)  XI 1A.Y1 1A.X12A.Y12A 

WRITE  (41,530)  Y11B.Y12B 

WRITE  (41.540)  X11C,Y11C,X12C.Y12C 

WRITE  (41,550)  XI 1G.Y1 1G.X12G.Y12C 

WRITE  (41,560)  Y11H.Y12H 

WRITE  (41,570)  X1112M,Y1112M,Z1112M 

PART  II  OF  THE  PROGRAM  CALCULATES  MOBILITY  FUNCTIONS 


Y11A  -  6.00*Y11A 
Y12A  -  6 . D0»Y12A 
Y11B  -  4.D0*Y11B 
Y12B  -  4.00*Y12B 
Y11C  -  8.00*Y11C 
Y12C  -  8.D0.Y12C 

AY(1 ,1)  -  Y11A 
AY(1 ,2)  -  Y12A 
AY(2, 1 )  -  Y12A 
AY(2,2)  -  Y11A 

AY(3.1)  -  Y11B 
AY(3,2)  -  Y12B 
AY(4.1)  -  -Y12B 
AY(4,2)  -  -Y11S 

AY(1 ,3)  -  Y11B 
AY(1 ,4)  -  -Y12B 
AY(2,3)  -  Y12B 
AY(2.4)  -  -Y11B 

AY(3,3)  -  Y11C 
AY(3,4)  -  Y12C 
AY(4,3)  -  Y12C 
AY(4,4)  -  Y11C 

DO  700  K-1 .4 
700  DU»rMY(K.1)  -  0.00 

CALL  DMAT1N(AY.4,0UM4Y,1,DETERM,4) 

Y11AM  -  AY(1 ,1)*6.D0 
Y12AM  -  AY(1 ,2)*6.00 
Y1 IBM  -  AY(3, 1 )*4.D0 
Y126M  -  AY(3,2)*4.00 
Y11CM  -  AY(3,3)*8.D0 
Y12CM  -  AY(3,4)*8.D0 


g 


ft 


t: 


i>: 


ft 


* 


’V. 

>; 

v. 

V. 


3 

3 


50 


XI 1AM 
XI 2AM 
X11CM 
X12CM 


X11A/(X11A*X11A  -  X12A»X12A) 
— X12A/(X1 1A»X1 1A  -  X12A*X12A) 
X11C/(X11C*X11C  -  X12C.X12C) 
-X12C/(X11C»X11C  -  X12C*X12C) 


Y11HM 

Y12HM 


-Y11C»Y11BM  -  Y12G»Y12BM  +  Y11H»Y11CM  +  Y12H»Y12CM 
Y1 1G*Y12BM  +  Y12G*Y1 IBM  4  Y11H*Y12CM  4  Y12H*Y1 1CM 


X11CM 

X12GM 

Y11GM 

Y12CM 


(X11G*X1 1AM  4  X12G*X12AM)/3.D0 
(XI 1G*X12AM  4  X12G*X11AM)/3.D0 

(Y11G*Y1 1AM  4  Y12G»Y12AM)/3.D0  -Y11H*Y11BM+Y12H*Y12BM 
(Y1 1G*Y12AM  4  Y12G*Y11AM)/3.D0  — Y1 1H*Y12BM+Y12H»Y1 IBM 


■X1112M  -  .  8D0* (XI 1G— X12G)» (X12GM-X1 1GM) 
•Y1112M  -  2 . 4D0» ( Y1 1G-Y1 2C) • (Y1 2GM-Y1 1GM) 
4  2.4D0«(Y11H+Y12H)«(Y11HMfY12HM) 

Z1WW  -  -Z1112M 


X1MM 

Y1M4 


C 

C 

C 

C 

C 

C 


END  OF  CALCULATIONS  FOR  THE  MOBILITY  FUNCTIONS 


COMPUTE  BATCHELOR  AND  GREEN'S  K.  L.  M  AND  J  FUNCTIONS 


».* 

BGK  - 

BGL  - 

BGM  -  - 

BGJ  -  - 

c 

1 

c 

WRITE  ( 

WRITE  ( 

WRITE  ( 

WRITE  ( 

WRITE  ( 

i 

WRITE  ( 
WRITE  ( 

5 

c 

WRITE  ( 

S00 

FORMAT 

> 

.■v 

510 

FORMAT 

S20 

FORMAT 

S30 

FORMAT 

5 

540 

FORMAT 

5S0 

FORMAT 

% 

o 

560 

FORMAT 

■s." 

570 

FORMAT 

0 

V 

620 

FORMAT 

& 

630 

FORMAT 

640 

FORMAT 

650 

FORMAT 

660 

FORMAT 

S 

670 

FORMAT 

sj 

680 

FORMAT 

Si 

s 

STOP 

END 

1 .00 


X11AM.Y11AM.X12AM.Y12AM 

Y11BM.Y12BM 

XI 1CM.Y1 1CM.X12CM.Y12CM 
X11GM.Y11GM.X12GM.Y12GM 
Y1 1W4.Y12HM 
XIKM.YIkM.ZIbM 
BGK.BGL.BGM.BGJ 


(15X. 'NUMBER  OF  POINTS:  '.14. 10X, 'R  -  ' .F8.3,//) 
(23X,'Xir .IIX.'YII'.IIX.'XIZ' .IIX.'YIZ*,/) 

(5X, 'RES.  FNCS.  A* ,F11 . 6.3F14.6,/) 

(5X, 'RES.  FNCS.  B’ .F25.6.F28.6,/) 

(SX.'RES.  FNCS.  C.F11.6.3F14.0./) 

(SX.’RES.  FNCS.  G' ,F1 1 .6.3F14.6,/) 

(5X, 'RES.  FNCS.  H* .F25.6.F2S.6./) 

(5X.’RES.  FCNS.  X11M+Xl2M,Y11M+Y12M,Z11MfZ12M: ' , 
F11.6.2F10.6.//) 

(5X. 'MOB.  FNCS.  A* ,F11 . 6.3F14.6,/) 

FNCS.  8' .F25.6.F28.6,/) 

FNCS.  C' , F11 . 6.3F14.6,/) 

FNCS.  G' ,F11 .6.3F14.6,/) 

FNCS.  H ' , F25 . 6 , F28 .  6 ,/) 

FNCS.  X1M,  YtM  ANO  Z1M: ' .F13.6.2F12.6,//) 
FCNS.  K.L.M  AND  J:\4F12.5) 


(SX.’MOB 
(SX.'MOB. 
(SX.’MOB. 
(SX.’MOB. 
(SX, 'MOB. 
(5X, 'MOB. 


SUBROUTINE  STOKES(R.NSTAR.M,SGN,NPT.RHS,AMN,BMN.CMN) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

dimension  A(see.see) .rhs(i } . ipvt(sm) ,z(see) . 

•  AMN(1).BhM(1),CMN(1) 

DATA  LDA/300/ 

C  DIMENSION  OF  MATRIX  A  MUST  EXCEED  OR  EQUAL  3*NPT. 

DM  -  DFLOAT (M) 

DO  10  K  -  1.NPT 

K2  -  K  +  NPT 

K3  -  K2  +  NPT 

N  -  K  4  NSTAR  -  1 

FN1  -  DFLOAT  (N+1 ) /DFLOAT  (4»W-2) 

FN2  -  DFLOAT(N-2)/DFLOAT(  N*(4»N-2)  ) 

DO  10  I  -  1.NPT 

12  -  I  4  NPT 

13  -  I  4  2*NPT 

THETA  -  DFLOAT ( 1—1 )/DFLOAT (NPT-1 ) *3. 141 S026S3S8979D0 
COI  -  DCOS(THETA) 

SI1  -  DSQRT ( 1 . D0-CO1 *C01 ) 

R2  -  DSQRT(SI1*SI1  4  (R4C01)**2) 

C02  -  — (R4C01  )/R2 
PNS1  -  PNS(N.M.COI) 

PNS2  -  PNS(N,M,C02) 

PPNS1  -  PNS(N.Mfl.COI) 

PPNS2  -  PNS(N.Mf1,C02) 

PNSP1  -  PNSP(N.M.COI) 

PNSP2  -  PNSP(N,M,C02) 

C  RATIOS  OF  SIN1/S1N2  FROM  FACTORING  OF  SINES  IN  EQUATIONS 
RZ  -  (  1.D0/R2  )**M 
RR  -  RZ/R2 
RPHI  -  RZ.R2 

A(I.K)  -  FN1 •(  C01*PNS1  -  RZ*SGN*C02*PNS2/R2«»N  ) 

O  -  FN2* (  PNSP1  -  RZ*SGN*PNSP2/R2**N  ) 

A(I ,K2)  -  -0FL0AT(N-M4l )•(  PNS(N4l .M.COI ) 

O  -  RZ«SGN*PNS(N4l,M.C02)/rt2«*(N42)  ) 

A(I ,K3)  -  DM»(  PNS1  -  R2*SGN*PNS2/R2**(N41)  ) 

C 

A(12,K)  -  (  FN1  4  FN2*DM  )•(  PNS1  4  RR*SGN*PNS2/R2«*N  ) 

O  4  FN2» (  C0l*PPNS1  4  RR»SGN*C02»PPNS2/R2*»N  ) 

A(I2.K2)  -  -OFLOAT(N+M4l)*(  PNS1  4  RR*SGN*PNS2/R2»«(N42)  ) 

•  -COI »PPNS1  -  RR*SCN*C02*PPNS2/R2**(N42) 

A(I2,K3)  -  -PPNS1  -  RR«SGN«PPNS2/R2»*(N4l) 

c 

A(I3.K)  -  — DM»FN2*(  PNS1  4  SGN*RPHI*PNS2/R2**N  ) 

A( I3.K2)  -  DM* (  PNS1  4  RPHI *SGN* PNS2/R2 • • ( N42 )  ) 

A(I3.K3)  -  PNSP1  4  RPHI*SGN»PNSP2/R2**(N41) 

10  CONTINUE 

C 

1DIM  -  3*NPT 

IF  (  M  .EQ.  0  )  IDIM  •  2«NPT 
CALL  DGECO(A, LDA, IDIM, IPVT,RCON,Z) 

CALL  DCESL(A, LDA. IDIM, IPVT.RHS.0) 

C 

DO  340  J-1.NPT 
J1  -  J  4  NPT 
J2  -  J1  4  NPT 
AMN(J)  -  RHS(J) 

BMN(J)  -  RHS(JI) 

340  CMN(J)  -  RHS(J2) 

RETURN 


SUBROUTINE  STOKESO(R.SGN.NPT.RHS.CWN) 

IMPLICIT  DOUBLE  PREC.-ION  (A-H.O-Z) 
dimension  A(iee.iee),RHsCi).CMN(i).iPVT(iee) 

DATA  LDA/100/ 

C  DIMENSION  OF  MATRIX  A  MUST  EXCEED  OR  EQUAL  W>T. 

DO  10  K  -  1.NPT 

N  -  K 

FN1  -  DFL0AT(N+l)/DFL0AT(4»N-2) 

FN2  -  OFLOAT (N-2)/DFLOAT (  N*(4.N-2)  ) 

DO  10  1  -  1.NPT 

THETA  -  OFLOAT ( 1-1 ) /OFLOAT (NPT-1 ) *3 . 1 4 1 59265358979D0 
COI  -  DCOS( THETA) 

SI1  -  DSORT ( 1 . D0-CO1 *C01 ) 

R2  -  DSQRT(SI1*SI1  4  (R4C01)**2) 

C02  -  -(R+C01)/R2 
PPNS1  -  PNS(N,1,C01) 

PPNS2  -  PNS(N,1 ,C02) 

C  RATIOS  OF  SIN1/SIN2  FROM  FACTOR INC  OF  SINES  IN  EQUATIONS 
RR  -  1.D0/R2 

A(I.K)  -  -PPNS1  -  RR.SGN.PPNS2/R2**(N+1) 

C 

10  CONTINUE 
C 

IDIM  -  NPT 

CALL  DGEFA(A.LDA. IDIM. IPVT, INFO) 

CALL  DGESL(A.LDA. IDIM. IPVT. RHS.0) 

DO  340  J-1.NPT 
340  CMN(J)  -  RHS(J) 

RETURN 

EM) 


DOUBLE  PRECISION  FUNCTION  PNS(N.M.X) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

IF  (M  .LE.  N)  CO  TO  1 
PNS  -  B.DB 
RETURN 
N1  ■  N-M 

IF  (M  .EO.  0)  CO  TO  IBB 
DM  -  DFLOAT(M) 

COEFF  -  1 .DB 
DO  5  1-1 ,M 

COEFF  -  COEFF*DFLOAT (2* 1—1 ) 

PB  -  COEFF 

IF  (N  .CT.  M)  CO  TO  7 

PNS  -  PB 

RETURN 

PI  .  PB*X*DFLOAT (2*M+1 ) 

IF  (N  .CT.  M+1)  CO  TO  B 
PNS  -  PI 
RETURN 
P  -  PI 
PM1NUS  -  PB 
DO  IB  1-2. N1 

PPLUS  -  2.DB*X*P  -  PMINUS 

+  DFLOAT(2*M>1)*(  X*P-PM1NUS  )/DFL0AT(I) 

PMINUS  -  P 
P  -  PPLUS 
PNS  -  PPLUS 
RETURN 

PB  -  1 .DB 

IF  (N  .CT.  M)  CO  TO  1B7 
PNS  -  PB 
RETURN 
PI  -  X 

IF  (N  .CT.  M+1)  CO  TO  IBB 

PNS  -  PI 

RETURN 

PMINUS  -  PB 

P  -  PI 

DO  110  1-2, N1 

PPLUS  -  2.DB*X*P  -  PMINUS  -  (X*P  -  PMINUS)/DFLOAT(I) 
PMINUS  -  P 
P  -  PPLUS 
PNS  -  PPLUS 
RETURN 
END 


iSSV/  .v.%  , 


DOUBLE  PRECISION  FUNCTION  PNSP(N.M.X) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-2) 

N1  -  N-M 

if  (m  .eo.  e)  co  to  iee 

DM  -  DFLOAT(M) 

COEFF  -  1 .DC 
DO  5  1-1 .M 

COEFF  -  COEFF »DFLOAT (2« 1-1 ) 

P9  -  COEFF 

PPB  -  -0M*X*P9 

IF  (N  .GT.  M)  CO  TO  7 

pnsp  «  ppb 

RETURN 

PI  -  P9»X»DFLOAT (2«M+1 ) 

PP1  -  DFLOAT (2*M+1 ) • (  P9»(1 .D9-X*»2)  +  X«PP9  ) 

IF  (N  .CT.  U+1)  GO  TO  9 

PNSP  -  PP1 

RETURN 

P  -  PI 

PP  -  PP1 

PMIN  -  P® 

DO  19  I-2.N1 

PPLUS  -  2.D9*X«P  -  PMIN  +  DFLOAT(2*U-1)*(X«P-PMIN)/DFLQAT(I) 
PPPLS  -  (  OFLOAT(I+2*M).(1.D9-X*«2)  +  DM*X*»2  )*P 
-  DM»X»PPLUS  +  X»PP 
PMIN  -  P 
P  -  PPLUS 
PP  -  PPPLS 
PNSP  -  PPPLS 
RETURN 

P6  -  1  .DO 

pp«  -  e.oe 

IF  (N  .CT.  •)  CO  TO  197 
PNSP  -  9.D9 
RETURN 
PI  -  X 

PP1  -  1 .09  -  X*X 

IF  (N  .CT.  1)  CO  TO  IBB 

PNSP  -  PP1 

RETURN 

PMIN  -  PB 

P  -  PI 

PP  -  PP1 

DO  119  I-2.N1 

PPLUS  -  2.DB«X«P  -  PMIN  -  (  X*P-PMIN  )/DFLOAT(I) 

PPPLS  -  DFL0AT(I)»PP1*P  ♦  X»PP 
PMIN  -  P 
P  -  PPLUS 
PP  -  PPPLS 
PNSP  -  PPPLS 
RETURN 
END 


.VA.VV.V.V 


SUBROUTINE  DMATIN(A,N,B,M,DETERM,NMAX) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  INDEX 

DIMENSION  A(NMAX,N) ,B(NMAX,M) ,PIV0T(166) , INDEX(IBB) 
DETERM-1. DB 
DO  26  1-1 .N 
PIVOT(I)-B.DB 
2B  INDEX(I)— B.DB 

DO  558  1-1 .N 
AMAX— B.DB 
DO  IBS  J-1.N 

IF  (PIVOT(J) .NE.B.OB)  CO  TO  1B5 
DO  IBB  K-1.N 

IF  (PIVOT(K) .NE.6.D0)  CO  TO  IBB 
TEMP-OABS(A(J.K)) 

IF  (TEMP. LT. AMAX)  CO  TO  IBB 
IROW-J 
ICOLUM-K 
AMAX-TEMP 
166  CONTINUE 
IBS  CONTINUE 

I NDEX ( I  )«4 . 896D3-DFL0AT ( I ROW)+OFLOAT ( ICOLUM) 

J-IROW 

AMAX— A(J . ICOLUM) 

C  SUPPRESS  CALCULATION  OF  DETERMINANT  •—••••••• 

C  DETERM-AMAX-DETERM 

IF  (DETERM. EQ.8.D6)  GO  TO  6BB 

PIVOT ( ICOLUM)— AMAX 

IF  (IRON. EQ. ICOLUM)  GO  TO  26B 

DETERM— OETERM 

DO  2BB  K— 1 ,N 

SWAP— A(J ,K) 

A(J. K)-A( ICOLUM, K) 

2BB  A( ICOLUM. K)-SWAP 

IF  (M.LE.B)  GO  TO  26B 
DO  2SB  K— 1 ,M 
SWAP— B(J ,K) 

B(J ,K)— B( ICOLUM, K) 

25B  B( ICOLUM.K)— SWAP 

266  K— ICOLUM 

A(ICOLUM.K)— 1 .DB 
DO  35B  K-1  ,N 

356  A( ICOLUM, K)—A( ICOLUM, K)/AMAX 

IF  (M.LE.B)  CO  TO  386 
DO  376  K-1 ,M 

376  B( ICOLUM, K)-B( ICOLUM, K)/AMAX 

386  DO  5S6  J-1 ,N 

IF  (J.EQ. ICOLUM)  CO  TO  556 
T-A(J, ICOLUM) 

A(J. ICOLUM) -B.DB 
DO  458  K-1 ,N 

456  A(J ,K)-A(J ,K)-A( ICOLUM, K)*T 

IF  (M.LE.B)  GO  TO  55B 
DO  5BB  K-1 ,M 

5BB  B(J ,K)-B(J , K)-B( ICOLUM, K)*T 

556  CONTINUE 
6BB  DO  716  1-1 ,N 
I1-N+1— I 

K-IO I NT( I NDEX ( I 1 )/4 . B96D3 ) 
IC0LUM-1D1NT(INDEX(I1)-4.696D3*DFL0AT(K)) 

IF  (K.EQ. ICOLUM)  GO  TO  716 
DO  765  J-1.N 
SWAP-A(J.K) 

A(J,K)-A(J, ICOLUM) 

765  A(J,ICOLUM)-SWAP 
716  CONTINUE 
RETURN 

EN°  V/. vV/v'. 


NUMBER  OF  POINTS:  12  R  -  4.000 


XII 

Y11 

X12 

Y12 

RES. 

FNCS. 

A 

1.169470 

1.043303 

-0.427212 

-0.204477 

RES. 

FNCS . 

B 

-0.020331 

0.098152 

|  RES. 

FNCS. 

C 

1.000295 

1 . 004226 

-0.015631 

0.008458 

|  RES. 

FNCS. 

C 

0.110028 

0.003564 

-0.255392 

-0.008486 

t  RES. 

FNCS. 

H 

-0.000331 

0.019701 

RES.  FCNS.  XI 1M+X12M.Y1 1M+Y12M.Z1 1M+Z12M:  1.094S51  0.971180  0.998097 


XII 

Y11 

X12 

Y12 

MOB.  FNCS. 

A 

0.986769 

0.999716 

0.360471 

0.195314 

MOB.  FNCS. 

B 

0.000120 

-0.031251 

MOB.  FNCS. 

C 

0.999948 

0.998914 

0.015625 

-0.007769 

MOB.  FNCS. 

G 

0.005504 

0.000020 

-0.070784 

-0.002604 

MOB.  FNCS. 

H 

-0.000750 

0.019569 

MOB.  FNCS. 

X1M. 

Y1M  AND  Z1M: 

-1.072550 

-0.970229 

-0.998097 

MOB.  FCNS.  K.L.M  AND  J:  -0.00190  -0.02787  0.16741  0.00184 


Sample  output  from  MANDR.FOR 
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